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ABSTRACT 



This thesis describes a Matlab computer program that implements Amott’s formulation of 
thermoacoustics [W.P. Amott, et. al., "General formulation of thermoacoustics for stacks 
having arbitrarily shaped pore cross sections," J. Acoustic. Soc. Am., 90(6), 3228-3237 
(1991)]. The program calculates the resonance frequency and the quality factor of a 
thermoacoustic prime mover below onset of self-oscilation. The results of this analysis are 
compared to measured values for both a closed end and an open end prime mover and to 
predictions of a standing wave analysis of prime movers [A. A. Atchley, "Standing wave 
analysis of a thermoacoustic prime mover below onset of self-oscillation," J. Acoustic. Soc. 
Am. ,92(5), 2907-2914 (1992)]. 
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I . INTRODUCTION 



Thermoacoustic heat transport is a process through which 
an acoustic field generates, or inversely is generated from, 
a flow of heat. Like every engine, a thermoacoustic engine can 
be configured as either type of classical heat engine, a heat 
pump (or refrigerator) or a prime mover. The purpose of this 
thesis is the computer implementation of a general formulation 
of thermoacoustics developed by Arnott et al [Ref. 1] using 
Matlab. In particular, we find the theoretical values of the 
quality factor and the resonance frequency of a prime mover as 
a function of the temperature difference applied across the 
prime mover stack. The program determines the impedance and 
normalized pressure distribution along the prime mover 
beginning from a known value of specific acoustic impedance 
and normalized pressure amplitude at one end. It finds the 
resonance frequency and quality factor through calculation of 
the pressure amplitude at the opposite end as a function of 
frequency. The result of the computer implementation are 
compared with measurements made with open and closed ended 
prime mover [Ref. 2,3], and with the results of a standing 
wave analysis of prime movers [Ref. 4]. 
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II. THEORY 



The goal of this chapter is to describe the basic 
geometry of a prime mover and to summarize the important 
points of Arnott's method. Parallel with this we explain how 
this theory is implemented in the computer program. The reader 
is referred to [Ref. 1] for full details of Arnott's method. 

A. ARNOTT'S METHOD. 

The necessary background for this thesis can be explained 
with the help of Fig.l which shows the basic construction of 
the thermoacoustic prime mover. It constists of a circular 
cross section acoustic resonator , rigidly capped at both 
ends. Situated within the resonator are the thermoacoustic 
elements consisting of an ambient heat exchanger, a prime 
mover stack and a hot heat exchanger. In our prime mover both 
the stack and heat exchangers are parallel plates spaced by a 
few thermal penetration depths. Complete details of the prime 
mover construction are given in [Ref. 2,3,4]. The traverse 
coordinates of the prime mover are taken to be x and y, and 
the longitudinal coordinate is z. The hot and cold sections 
are circular ducts open at the heat exchanger end and open or 
closed at the other end. The prime mover is divided into 
either 5 or 7 sections depending on the boundary conditions 
[Fig. 1,2]. The stack is divided into 11 subsections. 
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Figure 1. Rigid end Prime Mover 
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Figure 2. Open end Prime Mover 
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The computer program finds the normalized acoustic 
pressure and the specific acoustic impedance in all parts of 
the prime mover as a function of frequency. The pressure as a 
function of frequency gives the frequency response. From the 
frequency response we can determine the quality factor (Q) and 
the resonance frequency (f 0 ) . To do that we assume a linear 
temperature gradient across the stack via the hot and ambient 
(cold) heat exchangers. The stack is divided into eleven 
subsections. A constant temperature is assumed in each 
subsection of the stack. In all others parts of the prime 
mover we assume that the temperature is constant and equal to 
that of the nearest heat exchanger. The ambient temperature 
(T 0 ) is taken to be a function of z (T 0 (z)) only in the stack 
area . 

The program starts with specific acoustic impedance at 
the right end. There are two boundary conditions of interest: 
rigid and open. The specific acoustic impedance of a rigid end 
is : 



Z=(l+i) ^ 1 

2 (Y -1) 



P o c ' 



\ 



P pC d 

a) n 






(i) 



as described in [Ref. 5:p.529]. In the above equation N pr is 
the Prandtl number, c p is the isobaric heat capacity per unit 
mass, Y is the ratio c p /c v where c v is the constant volume heat 
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capacity per unit mass, T| is the viscosity and p 0 represents 
the ambient value of the density. The specific acoustic 
impedance of an open end is: 



Z=cp 0 ka ( 



ka 

4 



-0.6i) , 



( 2 ) 



where a is the radius of the tube and k=0)/c is the wavenumber, 
0) is the angular frequency and i=(-l) H [Ref. 6:p.202]. 

The acoustic pressure is specified at the same end. 
Because the Q is independent of the absolute value of the 
acoustic pressure in linear theory, the specified pressure is 
only a relative or normalized pressure. Next we find the 
pressure and the impedance at the entrance to the hot heat 
exchanger using the Rayleigh's impedance translation 
theorem [Ref. 1,5] . This is one difference between Arnott's and 
other methods. That is, other methods find the fluid variables 
by integrating the governing equations along the prime mover. 
The impedance translation theorem state that: within any 
homogeneous layer, with intrinsic (or characteristic) specific 
acoustic impedance Z lnt , the local specific impedance Z(z-l) at 
z-1 is related to that at z by [Ref. 5 :p. 139] : 



Z (y) cos (lei) -iZ lnC sin(JcI) 
int Z lnt cos (Icl) -iz(y) sin (icl) 
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Likewise, the acoustic pressure at the entrance to the hot 
heat exchanger is found through the pressure translation 
theorem: 



P 1 (z-1) =P 1 (z) {cos (kl) -i [ V/ ifl r 3 sin (Jcl) }, (4) 

Z(z) 

where Z(z) is the specific impedance at z, PjU) is the first 
order value of pressure, i.e, the acoustic pressure, and Z lnt 
is : 



Z - (5) 

int [Q F (A )k] ' 

In the above equation, Q represents the porosity and is equal 
to NA, where A is the cross sectional area of single pore and 
N is the number of pores per unit area, k is the complex wave 
number in the pore and given by: 

(k(A ,A T ) ) 2 =“i_ L_ ( Y -(Y -1)F(A r ) ) . (6) 

c 2 F(A ) 



The F (X) and F (X T ) factors arise from the equation for 
the z component of the acoustic velocity: 



u 



z (x,y, z) 



_ F(x,y,-A ) dP x {z) 
id) p 0 dz 



(7) 
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The F(x,y;X) satisfies the following differential equation 
[Ref. 1] : 

F(x,y;A ) +(-^— ) V^F(x,y;X )=1, (8) 

ik 2 

where the Laplacian operator V 2 T is given by: 



In our case where the pores of the stack are parallel plates, 
the F depends only on y and the shear wave number X or the 
dimensionless thermal disturbance number X T which are given 
by : 




( 10 ) 



and 



i 




( 11 ) 



The F(y;X) is given by: 



F(y,k )=1- 



cosh (/^T— ¥ ) 
2 a 



( 12 ) 



COSh( v / r l-y ) 
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In the above equation a is the half separation distance 
between the parallel plates of the stack or heat exchanger. 
The pore average of F(y;A) for heat exchangers and stack, is 
given by: 

F(A ) tanh(vFI^-) . (13) 

A 2 



Finally we use the boundary layer approximation in the 
circular regions of the prime mover: 

lim x _„F (A ) =1 - (5 /R) ( 1 +i ) , (14) 



In the above relation R is twice the hydraulic radius of the 
pore, or twice the ratio of the transverse pore area to the 
pore perimeter. Actually for heat exchangers and stack with 
parallel plates, R is the separation distance (2a) between the 
plates and for the circular parts it is the radius of the 
tube. 5 is the viscous or thermal penetration depth. 

Once Pj and Z are known at the entrance to heat exchanger, 
we need the values just inside the heat exchanger. Using the 
pressure and the impedance translation theorem we have a 
solution up to the stack where a linear temperature gradients 
exists. In the heat exchanger three things are different from 
the circular parts. These are: the wavenumbers in the heat 
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exchangers are complex, the porosity now is NA, and F (X) as we 
mention above is given by equation (13) . 

The point now is how to find pressure and impedance in 
the stack. In this case Arnott's method uses the following 
system of the first order differential equations that can be 
solved using Range-Kutta methods: 



dZ(z) =iie(z)Z inc (z) (1- Z(z)2 )+2a (z)Z(z), 
dz Z lnc (z) 2 



(15) 



and 




(16) 



where a(z) is: 




F( A T ) 



(17) 



T oz is the temperature gradient given by: 



, dr 0 (z) 
02 dz 



(18) 



(3 is the thermal expansion coefficient: 
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( 19 ) 



-( 



P = 



±£_ 

0 T 



P o 



p 



In our program, this system of equations is solved with the 
Runge-Kutta-Fehlberg method for every subsection of the stack. 
In so doing, the impedance and the pressure distribution in 
the stack is obtained. Now using the pressure and the 
impedance translation theorem we find the pressure and the 
impedance in the rest of the prime mover. 

Knowing the pressure at the left end we can calculate the 
quality factor, the resonance frequency and the maximum 
pressure amplitude. This is done by minimizing the following 
equation [Ref. 8]: 



( \p^y\£o )2 

Of 

2+_l 

£ fn O 2 



-ampli t 2 . 



( 20 ) 



where Q is the quality factor, |P(1) I is the maximum 
amplitude of pressure at the left end, f 0 is the resonance 
frequency, amp lit is the calculated amplitude of pressure at 
the left end and is a function of the frequency f. In the 
above relation there are two known variables f and amplit and 
three unknowns (Q,f 0 ,P(l)) that we want to find. To find these 
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unknowns, the 3-parameter least squares technique is used with 
the help of simplex method. 

Up to this point, Q and f 0 have been determined only by 
considering the pressure amplitude at the left end. We have 
not worried about matching the acoustic impedance at the left 
end. Impedance matching allows us to refine the values of Q 
and f 0 through the complex eigenf requency . Using this 
technique, one assumes the angular frequency is complex and 
given by: 



u=2nf 0 (l-l). (21) 

Substituting the initial value of f 0 and Q determined from the 
frequency response, into this equation gives an initial guess 
at the complex eigenf requency for the system. This frequency 
is used to calculate the impedance throughout the prime mover 
as explained earlier. The value of impedance in the gas at the 
left end is compared to the value for the boundary condition 
(i.e. the impedance of a rigid end with thermal loss). The 
difference between the computed and actual impedance is used 
to adjust the eigenf requency . The impedance is computed with 
the new eigenf requency until both real and imaginary parts of 
the impedance matche at the left end. The final values of f 0 
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and Q are determined from the real and imaginary parts of the 
eigenf requency that best matche the impedance . 



B. IMPLEMENTATION OF ARNOTT' S METHOD 

The program uses functions that already exist in Matlab 
where possible. For the Runge-Kutta-Fehlberg 5th order method 
we modified the existing function in Matlab to integrate with 
initial values of the integration variable greater than the 
final. This is necessary because we put the origin at the left 
end (z=0) and we begin our calculations from the right end 
(z=L) . One reason that we use this particular Runge-Kutta 
method is its accuracy. The Runge-Kutta-Fehlberg computes two 
Runge-Kutta estimates for the value y n+1 but of different 
orders of error. The global error in this numerical method is 
0 (h 5 ) and the local error is 0(h 6 ) [Ref. 7], 

To do the minimization, we use the function fmins with 
the help of simplex numerical method. The fmin function with 
simplex method is used to match the impedances at the left end 
using the technique of least squares. 
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III. THE MATLAB PROGRAM 



The goal of this chapter is to provide an overview of the 
structure of the Matlab program that implements Arnott's 
method. Furthermore it describes all the functions, special 
tricks and limitations of this code. It also explains how to 
run the program. 

A. OUTLINE OF THE MATLAB PROGRAM 

The program is divided into three parts. In the first 
part the variables are specified for every section of the 
prime mover. They are: the type of section (open, heat 
exchanger, stack) , the number of subsections that exists 
inside the section, the temperature at the left and right ends 
of the section, twice the hydraulic radius of the section, and 
the porosity of the section. Also, specified are the frequency 
range and the number of intervals in this range. 

The second part is the main program. It finds the 
impedance and the pressure distribution in every section of 
the prime mover for each frequency within the assumed range. 
To accomplish this it begins at the right end of the tube 
where it calculates the impedance of the end, rigid or open, 
and assigns an arbitrary value to the pressure. Continuing, it 
calculates using the translation theorems, the pressure and 
the impedance in all other sections except in the stack. In 
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the stack, where there exists a temperature gradient, the 
program calls the function "trode45" (a modified version of 
the function of Matlab "ode45") to integrate from the right 
end to the left end of the stack. This function solves 
differential equations using the 5th order Runge-Kutta 
numerical method with Fehlberg coefficients. 

The last part of the Matlab program uses the calculated 
pressure at the left end from the previous part to find the 
maximum amplitude, the resonance frequency and the quality 
factor. These parameters are found from a least square fit to 
equation (20) . The minimization is performed with the matlab 
function "fmins" with "options 1,2,3,14". Finally, the Q and 
f 0 estimates are used as a starting points for the final 
calculation of the quality factor and the resonance frequency 
using the complex eigenf requency method. This refinement 
usually is a small correction (in the fourth decimal place) to 
initial guesses. Finally, the program plots 1/Q and f D versus 
temperature difference and saves the results in a data file. 



B. LIMITATIONS OF THE PROGRAM, SPECIAL POINTS. 

There are two programs, one for the open end prime mover 
and one for the rigid end tube. To run the open end progam, 
one needs to run "arffl.m" with the functions "error, m- 
erro3 . m-argo3f . m" . To run the rigid end program one needs to 
run "argo2.m" with the functions "error . m-erro3 . m-argo3 . m" . 
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The difference in the two progams is that one uses "argo3.m" 
and the other uses "argof3.m". Before running one program or 
the other, one needs to switch these functions in the 
"erro3.m" function. 

Futhermore, there are some factors that influence the 
operation of the program. One important step is the choice of 
the initial range of the frequencies. To begin the program one 
needs to know the approximate resonance frequency as well as 
which longitudinal acoustic mode is to be used. Initially the 
frequency range is f 0 ± 20Hz. After the program runs for the 
first value of AT the frequency range is automatically 
adjusted to be: 




Another important factor is of course the number of 
points between the maximum and minimum frequency that the 
program uses. Usually this number is 500 but sometimes when 
the values of quality factor gets large this may have to be 
doubled (1000 points) . In these cases, it may also be 
necessary to use steps of 5 degrees Kelvin instead of the 
usual 10 degrees Kelvin. 

Finally, another important point concerns the "counters" 
that are used by the program. The most important counter is 
the "flag". If somebody wants to run the program for a case 
where the quality factor increases with AT, one needs to put 
the initial value of the "flag" counter equal to zero. 
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Otherwise, it should be set equal to one. The purpose of this 
counter is to force the program to scan for negative values of 
the inverse quality factor at the proper time. 
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IV. VALIDATION AND RESULTS 



In this chapter we compare the results of the program to 
experimental results for two different prime movers [Ref. 2 
and 3] as well as with the results of the standing wave 
analysis by Atchley [Ref 4]. 

The first results are for a rigid end prime mover 
described in Ref. 2 and 4. Comparison of both 1/Q and 
resonance frequency are made. Figures 3 through 5 show 1/Q vs 
AT (in degrees K) for three different mean gas pressures. The 
gas is helium. In these figures, the open circles represent 
measured values, the *'s represent the results of the standing 
wave analysis and the solid line represents the result of the 
Matlab program. The overall agreement between the present 
analysis and the measurements is good, but not as good as for 
the standing wave analysis. The present analysis agrees best 
at low values of AT. The reasons for the increased discrepancy 
at higher values of AT are not understood. One would expect 
the present analysis to agree better with measurements than 
the standing wave analysis. This is expected because the 
present analysis makes fewer assumptions. In the particular, 
it does a better job of taking into account changes in the 
acoustic pressure amplitude in differents parts of the prime 
mover. The standing wave analysis ignores any changes in 
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pressure amplitude due to changes in cross sectional area of 
the prime mover. 

Figures 6 through 8 show results of resonance frequency 
(in Hz) vs AT for the same three mean pressures. The agreement 
between the present analysis and the measured frequencies is 
within 1% in all cases. The standing wave analysis also agrees 
with the measured values to about the same accuracy. However, 
the two methods show different depedences on AT. This is 
because 1) the standing wave analysis considers only the 
temperature dependence of the sound speed, whereas the present 
analysis includes the effects of dispersion due to the 
boundary layers, and 2) the standing wave analysis ignores the 
finite impedance of the rigid ends. Of these two differences 
the first is the more important. This results in different 
curvatures of the two theoretical lines. 

The second set of comparisons is made for an open ended 
prime mover used in Che's thesis research [Ref. 3]. The 
unusual property of this prime mover is that while 1/Q 
decreases with AT for the first and third modes (as is usual 
for a prime mover) , 1/Q increases with AT for the second mode. 
This behavior provides a good test for theoretical models. The 
results for 1/Q vs AT are shown in figures 9 through 11 for 
the first, second and third modes, respectively. It is seen 
that the present analysis and the standing wave analysis are 
in close agreement in the first and second modes. The 
agreement diminishes slightly at the third mode. The agreement 
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with the measured values is worse than before. However, this 
is most likely due to contamination of the helium with air and 
water vapor as discuseed in Che's thesis. Another interesting 
aspect of the analysis is the prediction of an extremum in 1/Q 
for the second and third modes. 
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Figure 3. Graph of 1/Q vs AT (K) for the closed end prime 
mover filled with helium at 170 kPa mean pressure. 
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Figure 4. Graph of 1/Q vs AT (K) for the closed end prime 
mover filled with helium at 376 kPa mean pressure. 
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Figure 5. Graph of 1/Q vs AT (K) for the closed end prime 
mover filled with helium at 500 kPa mean pressure. 
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Figure 6. Graph of f 0 (Hz) vs AT (K) for the closed end prime 
mover filled with helium at 170 kPa mean pressure. 
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Figure 7. Graph of f 0 (Hz) vs AT (K) for the closed end prime 
mover filled with helium at 376 kPa mean pressure. 
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Figure 8. Graph of f 0 (Hz) vs AT (K) for the closed end prime 
mover filled with helium at 500 kPa mean pressure. 
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Figure 



9 . Graph 1/Q vs AT (K) for the first mode 
ended prime mover filled with helium gas 



of the 
at 101 



open 
kPa . 
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Figure 10. Graph of 1/Q vs AT (K) for the second mode of the 
open ended prime mover filled with helium gas at 101 kPa. 
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Figure 11. Graph of 1/Q vs AT (K) for the third mode of the 
open ended prime mover filled with helium gas at 101 kPa. 
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Figure 12. Graph of f 0 (Hz) vs AT (K) for the first mode of 
the open ended prime mover filled with helium gas 
at 101 kPa. 
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Figure 13. Graph of f Q (Hz) vs AT (K) for second mode of the 
open ended prime mover filled with helium gas at 101 kPa. 




Figure 14. Graph of f c (Hz) vs AT (K) for the third mode of 
the open ended prime mover filled with helium gas 
at 101 kPa. 
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V. SUMMARY 



The results of this thesis can be summarized as follows. 
Arnott's method was implemented in a Matlab progam. The 
impedance and the pressure distribution were calculated along 
the prime mover. From the frequency response at the left end 
and complex eigenf requency, the quality factor and the 
resonance frequency was estimated for both cases rigid and 
open end prime movers. The results of the analysis were 
compared with those of a standing wave analysis and with 
experimental results. In the most of the cases the results 
agree well. A future investigation can be the examination of 
the output work of the prime mover. 
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APPENDIX A. THE MAT LAB PROGRAM FOR THE OPEN END PRIME MOVER 



% Program for prime mover. 

% PART 1 (CHANGE) 
clear; 
clg; 



coun=0; 


% counter 


pt=0 ; 


% counter 


qual=zeros (45, 1) ; 


% quality factor 


inqual=zeros (45, 1) ; 


% inverse quality factor 


frO=zeros (45, 1) ; 


% resonant frequency 


deltaT=zeros (45, 1) ; 


% temperature difference 


point=0; 


% indicator 


f lag=l; 

slope=sum (' negative' ) ; 


% indicator 


load open.dat 


% standing wave analysis data 



% Main Loop for evrey temperature difference (itrgh) 
for it rgh=2 93 : -5 : 158 
coun=l+coun 



numfre=1000; 


% number of frequencies 


% establish some often used 


constants 


s=5 ; 

j=sqrt (-1) ; 


% number of sections 
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npr=2/3* (1+eps) ; 

gamma=5/3* (1+eps) ; 

cp= (2.5*8.3143/ (4 . e-3) ) * (1+eps) / 

ksolid=. 16* (1+eps) ; 

sqri= ( (1 + j) /sqrt (2) ) * (1+eps) ; 



pl=ones (18, numfre) . * (10~ (- 


-8) ) ; % presure inside the tube 


freq=zeros (1, numfre) / % 


f requncies 


w=zeros (1, numfre) ; % 


angular frequencies 


teren=zeros (1, s) ; % 


end of sections 


numblay=zeros (1, s) ; % 


layers of sections 


lengsec=zeros (1, s) ; % 


lengths of sections 


tempR=zeros (1, s) ; % 


rigth temperature of sections 


tempL=zeros (1, s) ; % 


left temperature of sections 


poros=zeros (1, s) ; % 


porocity of stack 



ratio=zeros ( 1 , s) ; % ratio 



% SET VALUES 




f remin=864 ; 


% minimum frequency 


f remax=904 ; 


% maximum frequency 


termin=sum (abs (' free' ) ) ; 


% end of the last section 


ampres= (1 . Ole +5) * (1+eps) ; 


% ambient presure 


dramp= (1 . e-16) ; 


% driver pressure 


% section "1" 

terenl=' opentu' ; 





teren (1,1) =sum (abs (terenl ) ) ; 
numblay (1, 1) =1; 
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lengsec (1, 1) = (66. 25e-2) * ( 1+eps) 
tempR(l, 1) =293* (1+eps) ; 
tempL (1,1) =293* (1+eps) ; 
ratio (1, 1) = (1 . 917e-2) * (1+eps) ; 
poros (1, 1) =1* (1+eps) ; 
ares (1, : ) =pi*ratio (1, : ) * (1+eps) 
% section "2" 

teren2 =/ hexch' ; 

teren (1,2) =sum (abs (teren2) ) ; 

numblay (1, 2) =1; 

lengsec (1, 2) = ( . 82e-2) * (1+eps) ; 
tempR(l, 2) =293* (1+eps) ; 
tempL (1,2) =293* (1+eps) ; 
ratio (l,2)=(5.08e-4) * (1+eps) ; 
poros (1, 2) =. 666* (1+eps) / 

% section "3" 
teren3=' stack' ; 
teren (1,3) =sum (abs (teren 3) ) ; 
numblay (1, 3) =11; 
lengsec (1, 3) = (2 . 59e-2) * (1+eps) ; 
tempR (1,3) =itrgh* (1+eps) ; 
tempL (1,3) =293* (1+eps) ; 
ratio (1, 3) = (8 . 636e-4) * (1+eps) ; 
poros (1,3)=. 89473* (1+eps) ; 

% section ”4" 
teren4=' hexch' ; 
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teren (1,4) =sum (abs (teren4) ) ; 
numblay (1, 4) =1; 

lengsec (1, 4) = ( . 82e-2) * (1+eps) ; 
tempR (1,4) =itrgh* (1+eps) / 
tempL (1,4) =itrgh* (1+eps) ; 
ratio (1, 4) = (5. 08e-4) * (1+eps) ; 
poros (1, 4) =. 666* (1+eps) ; 

% section "5" 

teren5=' opentu' ; 

teren (1,5) =sum (abs (teren5) ) ; 

numblay (1, 5) =1; 

lengsec (1,5) = (.6852)* (1+eps) ; 

tempR(l, 5) =itrgh* (1+eps) ; 

tempL (1,5) =itrgh* (1+eps) / 

ratio (1, 5) = (1 . 917e-2) * (1+eps) ; 

poros (1, 5) =1* (1+eps) ; 

% PART 2 
% MAIN PROGRAM 
% FREQUENCY DISTRIBUTION 

numtot=sum (numblay) +1 ; 
pt=pt+l ; 
if pt~=l 

df=frO (coun-1) /abs (qual (coun-1) ) ; 
f remin=f rO (coun-1 ) -df ; 
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fremax=frO (coun-1) +df ; 



end 

i=l : numfre; 
if i-=l 

f req=f remin; 
else 

freq (1 , i) =fremin+ (f remax- f remin ) . * (i . / ( numfre- 1) ) ; 

end 

w=2*pi*f req; 

% GET THE SPECIFIC ACOUSTIC IMPENDANCE AND PRESSURE AT ALL 
POINTS. 

% START AT THE RIGHT AN MOVE AT THE LEFT, 
dens (1 , numtot) =ampresM . Oe-3 / (tempR (1, s) *8 . 3143) ; % density 
vise (1, numtot) =1 . 887e-5* (tempR (l,s) /273. 15) A (. 6567) ; % 

viscosity 

kgas (1, numtot) =visc (1, numtot) *cp/npr; % termal conductivity 
sspeed (1 , numtot ) =972 . 8*sqrt (tempR (1 , s) /273 . 15) ; % speed 

% isothermal 

typel=' free' ; 
type2=' rigid' ; 

if termin==sum (abs (type2) ) 

% IMPEDANCE OF RIGID TERMINATION. 

fac (numtot, 1 rnumfre) =sqrt ( (dens (1, numtot) . . . 

* sspeed (1, numtot) A 2) . / (w. *visc (1, numtot) ) ) ; 
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z (numtot, 1 :numfre) = (1+ j) . * (dens (1, numtot) . *sspeed (1, numtot) . 
. . . *fac (numtot, 1 :numfre) *sqrt (npr) ) / (sqrt (2) * (gamma-1) ) ; 
% IMPEDANCE OF FREE TERMINATION 
else 

z (numtot, 1 rnumfre) =free (dens (1, numtot) , vise (1, numtot) , w, rati 
o (1, s) . . . , sspeed (1, numtot) , gamma, npr) ; 
end 

% IMPEDANCE TRANSLATE FOR THE OPEN TUBE OR HEAT EXCHANGER, 
f ault=' stack' ; 
truel=' opentu' ; 
true2=' hexch' ; 
for k=s : -1 : 1 
jup=0; 

jup=numtot-sum ( jup+numblay ( 1 , k : s) ) ; 
if teren (1, k) ~=sum (abs (fault ) ) 
dens (1, jup) =ampres*4 . Oe-3 / (tempR(l,k) . *8. 3143) ; % density 
visc(l, jup)=1.887e-5.* (tempR (1 , k) /273.15) . * ( • 6567) ; 

% viscosity 

kgas ( 1 , jup) =visc ( 1 , jup) . *cp/npr; % termal conductivity 
sspeed ( 1 , jup) =972 . 8 . *sqrt (tempR ( 1 , k) /273 . 15) ; % speed 

% isothermal 

% GET LAMBDA AND LAMBDA (T) 

lambda (jup, 1 :numfre) =ratio (1, k) . *sqrt (dens (1, jup) . *w. /vise (1 
, jup) ) ; 

lambdt ( jup, 1 : numfre) =sqrt (npr) . *lambda ( jup, 1 : numfre) ; 
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% GET F (L) , F (L (T) ) AND WAVENUMBERS FOR OPEN TUBE 

PARTS 

if teren (1, k) ==sum(abs (truel) ) 
f lam ( jup, 1 : numf re) =1- (1+ j) . *sqrt (2) . /lambda ( jup, 1 : numfre) ; 

% f(l) 

f lamt (jup, 1 : numfre) =1- ( 1+ j ) . *sqrt (2) . /lambdt ( jup, 1 : numfre) ; 

% f (1 (t) ) 

facl = ( 1+ ( gamma- 1) /sqrt (npr) ) / sqrt (2) ; 
ka ( jup, 1 : numf re) = (w. /sspeed (1, jup) ) . . . 

. * (1+ (1+ j) . *facl. /lambda (jup, 1: numfre) ) ; 
else 

% GET F (L) , F (L (T) ) AND WAVENUMBERS FOR HEAT 

EXCHANGERS TYPE "SLIT". 

sqrmi= ( 1 — j ) /sqrt (2) ; 

ar ( jup, 1 : numfre) =sqrmi . *lambda (jup, 1 : numfre) ; 
argum ( jup, 1 : numfre) =exp (-ar ( jup, 1 : numfre) ) ; 
ar ( jup, 1 : numfre) =ar ( jup, 1 : numfre) /2; 
ctanh ( jup, 1 : numfre) = (1-argum ( jup, 1 : numfre) ) . / ( 1+argum (jup, 1 : 
numfre) ) ; 

flam ( jup, 1 : numfre) =l-ctanh ( jup, 1 : numfre) . /ar ( jup, 1 : numfre) ; 

% f(l) 

ar ( jup, 1 : numfre) =sqrmi*lambdt (jup, 1 : numfre) ; 
argum ( jup, 1 : numfre) =exp (-ar (jup, 1 : numfre) ) ; 
ar ( jup, 1 : numfre) =ar ( jup, 1 : numfre) / 2; 
ctanh ( jup, 1 : numfre) = (1-argum ( jup, 1 : numfre) ) . / (1+argum (jup, 1 : 
numfre) ) ; 
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f lamt (jup, 1 : numf re ) =l-ctanh ( jup, 1 : numf re) . /ar ( jup, 1 : numfre) ; 

% f (1 (t) ) 



low=w/ sspeed (1 , jup) ; 

ka (jup, 1 : numfre) =low . *sqrt ( (gamma- ( gamma- 1 ) . . . 

. * f lamt ( jup, 1 : numfre) ) . / flam (jup, 1 : numfre) ) ; 
end 

% TRANSLATION THEOREM 

comden ( jup, 1 : numfre) =dens (1 , jup) . /flam ( jup, 1 : numfre) ; 
zint ( jup, 1 : numfre) =comden ( jup, 1 : numfre) .*w./(ka(jup,l: numfre 
) . *poros (1, k) ) ; 

dsub (1, k) =lengsec (1, k) . /numblay (1, k) ; 
ar ( jup, 1 : numfre) =ka ( jup, 1 : numfre) . *dsub ( 1 , k) ; 
sn ( jup, 1 : numfre) =sin (ar ( jup, 1 : numfre) ) ; 
cs ( jup, 1 : numfre) =cos (ar ( jup, 1 : numfre) ) ; 
ct (jup, 1 : numfre) =cs ( jup, 1 : numfre) . /sn ( jup, 1 : numfre) ; 
fac3 ( jup, 1 : numfre) =zint ( jup, 1 : numfre) . /z ( jup+1 , 1 : numfre) ; 

z (jup, 1 : numfre) =zint ( jup, 1 : numfre) . * (ct ( jup, 1 : numfre) . . . 
- j . *fac3 ( jup, 1 : numfre) ) ,/(fac3(jup,l: numfre) . *ct ( jup, 1 : numf r 
e) - j) ; 

pi (jup, 1 :numfre) =pl (jup+1, 1 :numfre) . * (cs (jup, 1 rnumfre) . . . 
-j.*fac3(jup,l: numfre) . *sn ( jup, 1 : numfre) ) ; 
else 

% STACK WITH LINEAR DISTRIBUTION OF TEMPERATURE BETWEEN 
THE ENDS OF THE STACK 

toz= (tempR ( 1 , k) -tempL ( 1 , k) ) /lengsec (1 , k) ; % approximation 

dz/dt 
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ns=0; 



for la=k+numblay (1, k) -1 : -1 : k 
ns=ns+l ; 

tens=tempR ( 1 , k) - (tempR ( 1 , k) -tempL ( 1 , k) ) * ns /numb lay ( 1 , k) ; 
tnsml=tempR ( 1 , k) - (tempR ( 1, k) -tempL ( 1 , k) ) * (ns-1) / numb lay ( 1 , k) 



tave ( 1 , la) = (tens+tnsml ) /2; 

dsub (1, k) =lengsec (1, k) /numblay (1, k) ; 

dens (1, la) =ampres*4 . Oe-3 / (tave (1, la) *8. 3143) ; % density 

vise (1, la) =1 . 887e-5* (tave (1, la) /273 . 15) ^ ( • 6567) ; % viscosity 
kgas (1, la) =visc (1, jup) *cp/npr; % termal conductivity 
sspeed (1, la) =972 . 8*sqrt (tave (1, la) /273 . 15) ; % speed 

% isothermal 



% GET LAMBDA AND LAMBDA (T) 

lambda (la, 1 : numfre) =ratio (1, k) . *sqrt (dens (1,1a) .*w./visc(l,l 
a) ) ; 

lambdt (la, 1 : numfre) =sqrt (npr) *lambda (la, 1 : numfre) ; 
ar (la, 1 : numfre) =sqrmi*lambda (la, 1 : numfre) ; 
argum (la, 1 :numfre) =exp (-ar (la, 1 :numfre) ) ; 
ar (la, 1 : numfre) =ar (la, 1 : numfre) / 2; 
ctanh (la, 1 : numfre) = (1 -argum (la, 1 : numfre) ) . / (1+argum (la, 1 : num 
fre) ) ; 

flam (la, 1 : numfre) =1 -ctanh (la, 1 : numfre) . /ar (la, 1 : numfre) ; 

% f (1) 

ar ( la, 1 : numfre ) =sqrmi * lambdt (la, 1 : numfre) ; 
argum (la, 1 : numfre) =exp (-ar (la, 1 : numfre) ) ; 
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ar (la, 1 rnumfre) =ar (la, 1 rnumfre) / 2; 

ctanh (la, 1 : numf re) = (1-argum (la, 1 : numf re) ) . / (1+argum (la, 1 : num 
fre) ) / 

flamt (la, 1 : numf re) =1 -ctanh (la, 1 rnumfre) . /ar (la, 1 rnumfre) ; 

% f (1 (t) ) 

low=w/sspeed (1 , la) ; 

ka (la, 1 : numf re) =low. *sqrt ( (gamma- (gamma-1) . . . 

. * flamt (la, 1 : numf re) ) . /flam (la, 1: numf re) ) ; 
zint (la, 1 rnumfre) =dens (l,la).*w... 

. / (poros (1, k) .*flam(la,l rnumfre) . *ka (la, 1 rnumfre) ) ; 
alprim (la, 1 r numf re) =toz* (flamt (la, 1 : numf re) . . . 

. /flam (la, 1 rnumfre) -1) . / (2*tave (1, la) * (1-npr) ) ; 
kal=ka(la, 1 rnumfre) . ' ; 
zintl=zint (la, 1 : numf re) . ' ; 
alprl=alprim (la, 1 : numf re) . ' ; 

zetaO= (lengsec (1, k) - (ns-1) * lengsec (1, k) /numblay (1, k) ) ; 
zetaf in= ( lengsec ( 1 , k) - (ns) * lengsec ( 1 , k) /numblay ( 1 , k) ) ; 
zO=z (la+1, 1 r numf re ) . ' ; 
tol=l . e-4 ; 

[zeta, zout] =trode4 5 (' tube' , zetaO, zetaf in, zO, tol, kal, zintl, al 
prl) ; 

n=length (zeta) ; 

z (la, 1 r numf re) =zout (n, r ) ; 

zdumy=zO; 

plO=pl (la+1, 1 r numf re) . ' ; 
tol=l . e-4 ; 
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[ zeta, pi out ] =t rode45 (' tubel ' , zetaO, zetafin,plO,tol,kal,zintl 
, zdumy ) ; 

nl=length (zeta) ; 

pi (la, 1 :numfre) =plout (nl, : ) ; 

end 

end 

end 



% MINIMIZATION TO GET RESONANCE FREQUENCY AND QUALITY 
FACTOR. 

tfun=- j*w. *z (1, 1 :numfre) *dramp. /pi (1, 1 : numfre) ; 
for a=l : s 

pi (a, 1 : numf re) =pl (a, 1 : numf re) . *tfun; 
end 

amplit=zeros (1, numfre) ; 

amplit (1,1: numfre ) =abs (pi ( 1 , 1 : numfre) ) ; 

[maxamp, i ] =max (amplit ) ; 
resf re=f req (i ) 
amphaf=maxamp/sqrt (2) ; 
q=0; 

f rehaf=0; 
if point==0 

for s=2:numfre 

if amplit (s) >=amphaf & amplit (s-1 ) <=amphaf 
f rehaf= (freq (s) +f req (s-1 ) ) /2; 
q=0 . 5*resf re/ ( resf re-f rehaf ) 
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end 



end 

else 

for s=2:numfre 

if amplit (s) <=amphaf & amplit (s-1 ) >=amphaf 
ama=l 

f rehaf= (f req (s) +f req (s-1 ) ) / 2 
q=0 . 5*resf re/ (resf re-f rehaf ) 
end 

end 

end 

if q~=0 

xO=zeros (3,1); 
xO ( 1 , 1 ) =maxamp; 
xO (2, 1) =resfre; 
xO (3, 1) =q; 
options (1) =0; 
options (2) =1 .e-4; 
options (3) =1 .e-6; 
options (14) =2000; 

est=fmins ( ' erro' , xO, options, [ ] , amplit, f req) ; 

inqual (coun, 1) =est (3, 1) ~ (-1) ; 

frO (coun, 1) =est (2, 1) ; 

deltaT (coun, 1) =abs (itrgh-293) ; 

qual (coun, 1) =est (3, 1) ; 

end 
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if coun~=l 



dif=inqual (coun-1, 1) -inqual (coun, 1) ; 
if dif<0 & coun==2 

slope=sum (abs ( ' positive' ) ) ; 
elseif dif<0 & coun==3 

slope=sum (abs ( ' positive' ) ) ; 

end 

if flag==0 & slope~=sum (abs (' positive' ) ) 
if dif <=0 | q==0 
point=l ; 
f lag=l; 

for s=2:numfre 

if amplit (s) <=amphaf & amplit (s-1) >=amphaf 
f rehaf= (f req (s) +f req (s-1 ) ) 12 ; 
q=0 . 5*resf re/ (resf re-f rehaf ) 
end 
end 

xO ( 1 , 1 ) =maxamp; 
xO (2, 1 ) =resfre; 
xO (3, 1) =q; 
options ( 1 ) =0; 
options (2) =1 . e-4 ; 
options (3) =1 .e-6; 
options (14) =2000; 

est=fmins ( ' erro' , xO, options, [ ] , amplit, f req) ; 
inqual (coun, 1) =est (3, 1) A (-1) ; 
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frO (coun, 1) =est (2, 1) ; 
deltaT (coun, 1) =abs (itrgh-293) ; 
qual (coun, 1) =est (3, 1) ; 

end 

end 

end 

% PART THREE 

% MINIMIZATION USING COMPLEX ANGULAR FREQUENCIES, 
options (1) =0; 
options (2) =1 . e-4; 
options ( 3) =1 . e-4 ; 
options ( 14 ) =700; 

wl=2*pi*fr0 (coun, 1) * ( 1 — j / (2*qual (coun, 1) ) ) ; 
y2=wl+wl*10 . e-6; 
yl=wl-wl*10 .e-6; 

ella=fmin (' erro3' , yl, y2, options, itrgh) ; 
frO (coun, 1) =real (ella) / (2*pi) ; 

qual (coun, l)=-real(ella) / (2*imag (ella) ) ; 
inqual (coun, 1) =qual (coun, 1) A (-1) ; 
end 

% PLOTS AND STORE THE DATA (CHANGE) . 

y= ( (inqual (1 : coun, 1 ) ) ' ; (deltaT (1 : coun, 1) ) ' ; 

( f rO ( 1 : coun, 1) ) ' ] ; 

fid=fopen (' inqf3.txt' , ' w+' ) ; 
fprintf (f id, ' Inqual=%-9 . 6f T=%-4.1f RF=%-5 . 2 f \n' , y ) ; 
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load compf2.dat 
figure 



plot (deltaT (1 : coun) , inqual (1 : coun) , compf 2 ( : , 1 ) , compf2 ( : , 2) , ' 
o' . . . 

, open ( : , 1 ) , open ( : , 4 ) , ' * ' ) 
title (' Inverse Quality Factor V.S Temperature 
Difference . Open End'); 

grid 

xlabel (' Negative Delta "T" '); 
y label (' 1/Q # ) ; 

gtext (' Ambient pressure 101 Kpa' ) 

gtext (' Experimental Result "o”') 
gtext (' Computer Result 
gtext ('Third Mode') 
print -deps inquf3 
figure 

plot (deltaT (l:coun),fr0(l: coun) ) 
xlabel (' Negative Delta "T" '); 
ylabel ('Resonance Frequence' ) ; 
title (' Resonance Frequency V.S Temperature 
Difference . Open End'); 
gtext ('Ambient pressure 101 Kpa') 
gtext ('Third Mode') 
grid 

print -deps frf3 
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% FUNCTION TRODE45.M 



function . . . 



[tout, yout ] =trode45 (FunFcn, tO,tfinal,yO,tol,kal, zintl, alprl ) 
alpha= [1/4 3/8 12/13 1 1/2]'; 



beta= [ [1 0 000 0 ] / 4 

[3 9 0 0 0 0 ] / 3 2 

[ 1932 -7200 7296 0 0 OJ/2197 

[ 8341 -32832 29440 -845 0 0 ] / 4 1 0 4 

[-6080 41040 -28352 9295 -5643 0J/20520]'; 

gamma= [ [902880 0 3953664 3855735 -1371249 277020 ] /7618050 
[ -2090 0 22528 21970 -15048 -27360 ] /752400 

]'; 

% initialization 



t=t0; 

hmax= (t-tfinal) /5; 
hmin= (t-tfinal) /2000; 

h= (t-tfinal) /100; 
y=y0 ( : ) ; 



f=y*zeros (1,6); 
tout=t ; 
yout=y . ' ; 

tau=tol*max (norm (y, ' inf ' ) , 1) ; 

% Main loop 

while (t>tfinal) & (h>=hmin) 
if t-h>tfinal ,h=t-tfinal; end 
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% Compute the slopes 
temp=feval (FunFcn, kal, zintl, alprl, t , y) ; 
f ( : , 1 ) =temp ( : ) ; 
for j=l:5 

temp=feval (FunFcn, kal, zintl, alprl, t-alpha ( j) *h, y-h*f *beta ( 

j) > ; 

f ( : , j + 1) =temp ( : ) ; 
end 

% Estimate the error and the acceptaple error 
delta=norm (h*f *gamma ( : , 2) , ' inf' ) ; 
tau=tol*max (norm (y, ' inf' ) , 1 ) ; 

% Update the solution only if the error is acceptaple 
if delta<=tau 
t=t-h; 

y=y-h*f *gamma ( : , 1 ) ; 
tout= [tout ; t ] ; 
yout= [yout;y . ' ] ; 
end 

% Update the step size 
pow=l/5; 
if delta~=0.0 

h=min (hmax, 0 . 8*h* (tau/delta) A pow) ; 

end 

end 

if (tfinalct) 

disp( 'singularity likely.') 
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t 

end 

% FUNCTION "TUBE . M" GIVES THE DERIVATIVE OF IMPEDANCE 
function dzdz=tube (kal, zintl, alprl, zeta, z) 
dzdz=(j.*kal.*zintl.*(l~(z. /zintl) . A 2)+2.*alprl.*z) ; 

% FUNCTION "TUBEl.M" GIVES THE DERIVATIVE OF PRESSURE 
function dpdz=tubel (kal, zintl, zdumy, zeta, pi) 
dpdz= j . *pl . *kal . *zintl . /zdumy; 

% FUNCTION "ERRO.M" FOR "FMINS" . LEAST SQUARS TECHNIQUE. 

function rmserr=erro (xO, amplit, freq) 
mamp=xO (1,1); 
f 0=x0 (2,1) ; 
z=xO (3,1); 

focl= (mamp*fO) . / (z*freq) ; 

err= ( (foe 1.^2) ./ ((fO. /f req-f req/f 0 ) . / '2 + l/z / '2) ) -amplit . / '2; 
rmserr=sum (err . ~2) ; 

% FUNCTION "ERR03.M" FOR "FMIN" . LEAST SQUARS TECHNIQUE 
function fplus=erro3 (wl , itrgh) 

[ z , dens, s speed, vise] =argo3f (wl, itrgh) ; 

gamma=5/3* (1+eps) ; 

npr=2/3* (1+eps) ; 

sqri= ( (1+ j) / sqrt (2) ) * (1+eps) ; 
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f ofw=sqri *sqrt (dens A 3*sspeed A 4*npr . . . 

. / (wl . *visc) ) . / (gamma-1) ; 
sl=real ( fofw) ; 
s2=real (z) ; 
s3=imag (fofw) ; 
s4=imag (z) ; 

fplus = (sl A 2-s2 A 2) A 2 + (s3 A 2-s4 A 2) A 2; 



% FUNCTION "ARG03F.M" . "USED BY ERR03.M" 
function 

[z, dens, s speed, vise] =argo3 (wl, itrgh) 
w=wl ; 

% establish some often used constants 

s=5; % number of sections 

j=sqrt (-1) ; 
npr=2/3* (1+eps) ; 
gamma=5/3* (1+eps) ; 
cp= (2.5*8. 3143/ (4.e-3) ) * (1+eps) ; 
ksolid=. 16* (1+eps) ; 
sqri= ( (1 + j) / sqrt (2) ) * (1+eps) ; 
teren=zeros (s, 1) ; 
numblay=zeros (1, s) ; 
lengsec=zeros (1, s) ; 
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tempR=zeros (1, s) ; 
tempL=zeros ( 1 , s) ; 
poros = zeros (1, s) / 
ratio=zeros (1, s) ; 

% SET VALUES 

termin=sum (abs ( ' free' ) ) ; 
ampres= (1 . Ole + 5) * (1+eps) ; 
dramp= ( 1 . e-8) ; 



% section "1" 

terenl = ' opentu' ; 

teren (1,1) =sum (abs (terenl ) ) ; 

numblay ( 1 , 1 ) =1 / 

lengsec (l,l) = (65.34e-2)* (1+eps) 
terapRd, 1) =293* (1+eps) ; 
tempL (1, 1) =2 93* (1+eps) ; 
ratio (1, 1) = (1 . 917e-2) * (1+eps) ; 
poros (1,1) =1* (1+eps) ; 
ares (1, : ) =pi* ratio (1, : ) * (1+eps) 
% section "2" 
teren2=' hexch' ; 
teren (1,2) =sum (abs (teren2) ) / 
numblay ( 1 , 2) =1; 

lengsec (l,2) = (.82e-2) * (1+eps) ; 
tempR(l, 2) =293* (1+eps) ; 
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tempL ( 1 , 2 ) =293* (1+eps) ; 
ratio (1 , 2) = (5 . 08e-4 ) * (1+eps) ; 
poros (1, 2) =. 666* (1+eps) ; 

% section "3" 
teren3=' stack' ; 
teren (1,3) =sum (abs (teren3) ) ; 
numblay (1, 3) =11; 
lengsec (1,3) = (2.59e-2) * (1+eps) ; 
tempR(l, 3) =itrgh* (1+eps) ; 
tempL (1 , 3) =293* (1+eps) ; 
ratio ( 1, 3) = (8 . 636e-4 ) * (1+eps) ; 
poros (1,3)=. 89473* (1+eps) ; 

% section "4" 
teren4 = ' hexch' ; 
teren (1,4) =sum(abs (teren4) ) ; 
numblay (1, 4) =1; 

lengsec (l,4) = (.82e-2)* (1+eps) / 
tempR (1,4) =itrgh* (1+eps) ; 
tempL (1,4) =itrgh* (1+eps) ; 
ratio (1, 4) = (5. 08e-4) * (1+eps) ; 
poros (1, 4) =. 666* (1+eps) ; 

% section "5” 

teren5=' opentu' ; 

teren (1,5) =sum (abs (teren5) ) ; 

numblay ( 1 , 5) =1 ; 

lengsec (1,5) = (.6852) * (1+eps) ; 
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tempR (1, 5) =itrgh* (1+eps) ; 
tempL (1,5) =itrgh* (1+eps) ; 
ratio (1, 5) = (1 . 917e-2) * (1+eps) ; 
poros (1, 5) =1* (1+eps) ; 
if s>=7 
% section "6" 

teren6=' hexch' ; 

teren (1,6) =sum (abs (teren6) ) ; 

numblay (1, 6) = 1; 

lengsec (1, 6) = (1 . 02e-2) * (1+eps) ; 
tempR(l, 6) =293* (1+eps) ; 
tempL (1, 6) =2 93* (1+eps) ; 
ratio (1, 6) = (1 . 02e-3) * (1+eps) ; 
poros (1, 6) =. 667* (1+eps) ; 

% section "7" 

teren7=' opentu' ; 

teren (1,7) =sum (abs (teren7) ) ; 

numblay ( 1 , 7) =1 ; 

lengsec (1, 7) = (5 . 483e-2) * (1+eps) ; 
tempR (1, 7) =293* (1+eps) ; 
tempL (1, 7) =293* (1+eps) ; 
ratio (1, 7) =1 . 917e-2* (1+eps) ; 
poros (1, 7) =1* (1+eps) ; 
end 

% START AT THE RIGHT AND MOVE AT THE LEFT. 
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numtot=sum (numblay) +1; 

dens (1, numtot) =ampres*4 . Oe-3 / (tempR(l, s) *8 . 3143) ; % density 
vise (1, numtot) =1.887e-5* (tempR (1, s) /273.15) ~ (.6567) ; 

% viscosity 

kgas (1, numtot) =visc (1, numtot) *cp/npr; % termal conductivity 
sspeed ( 1, numtot ) =972 . 8*sqrt (tempR (1 , s) /273 . 15) ; % speed 

% isothermal 

typel=' free' ; 
type2=' rigid' / 

if termin==sum (abs (type2 ) ) 

% IMPEDANCE OF RIGID TERMINATION. 
fac3 (1, numtot) =sqrt ( (dens (1, numtot) . . . 

*sspeed (1, numtot) . A 2) . / (w. *visc (1, numtot) ) ) ; 
z (1, numtot) = (1+ j) . * (dens (1, numtot) . *sspeed (1, numtot) . . . 

. *fac3 (1, numtot) *sqrt (npr) ) / (sqrt (2) * (gamma-1) ) ; 

% IMPEDANCE OF FREE TERMINATION 
else 

z (1, numtot) =free (dens (1, numtot) , vise (1, numtot) , w, ratio (1, s) . 
, sspeed (1, numtot) , gamma, npr) ; 

end 

% IMPEDANCE TRANSLATE FOR THE OPEN TUBE OR HEAT EXCHANGER, 
f ault=' stack' ; 
truel=' opentu' ; 
true2=' hexch' ; 
for k=s:-l:l 
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jup=0 ; 

jup=numtot-sum ( jup+numblay ( 1 , k : s) ) ; 
if teren (1, k) ~=sum (abs (fault) ) 
dens (1, jup) =ampres*4 . Oe-3 / (tempR ( 1 , k) . * 8 . 3143) ; % density 
vise (1, jup) =1.887e-5.* (tempR (l,k)/273.15) . ~ ( • 6567) ; 

% viscosity 

kgas (1, jup) =visc (1, jup) . *cp/npr; % termal conductivity 
sspeed ( 1 , jup) =972 . 8 . *sqrt (tempR (1, k) /273 . 15) ; % speed 

% isothermal 



% GET LAMBDA AND LAMBDA (T) 

lambda (1, jup) =ratio (1, k) . *sqrt (dens (1, jup) . *w. /vise (1, jup) ) ; 
lambdt (1, jup) =sqrt (npr) . * lambda (1, jup) ; 

% GET F (L) , F (L (T) ) AND WAVENUMBERS FOR OPEN TUBE 

PARTS 

if teren ( 1 , k) ==sum (abs (truel ) ) 
flam (1, jup) =1- (1 + j) . *sqrt (2) . /lambda (1, jup) ; % f < 1 ) 

flamt (1, jup) =1- (1 + j) . *sqrt (2) . /lambdt (1, jup) ; % f(l(t)) 

fac31= (1+ (gamma-1) / sqrt (npr) ) / sqrt (2) ; 
ka ( 1, jup) = (w . /sspeed (1 , jup) ) . . . 

. * (1+ (1 + j) .*fac31. /lambda (1, jup) ) ; 

else 

% GET F (L) , F (L (T) ) AND WAVENUMBERS FOR HEAT 

EXCHANGERS TYPE "SLIT". 

sqrmi= (1-j) /sqrt (2) ; 

ar (1, jup) =sqrmi . *lambda (1, jup) ; 
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argum (1, jup) =exp (-ar (1, jup) ) ; 
ar (1, jup) =ar (1, jup) 12 ; 

ctanh (1 , jup) = ( 1-argum (1 , jup) ) . / ( 1+argum ( 1 , jup) ) ; 
flam(l, jup) =l-ctanh (1, jup) . /ar (1, jup) ; % f (1) 

ar (1, jup) =sqrmi*lambdt (1, jup) ; 
argum (1, jup) =exp (-ar (1, jup) ) ; 
ar (1, jup) =ar (1, jup) /2; 

ctanh (1, jup) = (1-argum (1, jup) ) . / (1+argum (1, jup) ) ; 
f lamt (1, jup) =l-ctanh (1, jup) . /ar (1, jup) ; % f(l(t)) 
low=w/sspeed ( 1, jup) / 

ka (1, jup) =low. *sqrt ( (gamma- (gamma-1) . . . 

. *flamt (1, jup) ) . /flam (1, jup) ) ; 

end 

% TRANSLATION THEOREM 
comden (1, jup) =dens (1, jup) . /flam (1, jup) ; 
zint (1, jup) =comden (1, jup) . *w. / (ka (1, jup) . *poros (1, k) ) ; 
dsub ( 1 , k) =lengsec ( 1 , k) . /numb lay (1 , k) ; 
ar (1, jup) =ka (1, jup) . *dsub (1, k) ; 
sn (1, jup) =sin (ar (1, jup) ) ; 
cs (1, jup) =cos (ar (1, jup) ) ; 
ct (1, jup) =cs (1, jup) . /sn (1, jup) ; 
fac3 (1, jup) =zint (1, jup) ./z (1, 1 + jup) ; 
z (1, jup) =zint (1, jup) . * (ct (1, jup) . . . 

- j. *fac3 (1, jup) ) . / (fac3 (1, jup) . *ct (1, jup) — j ) / 

else 
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% STACK WITH LINEAR DISTRIBUTION BETWEEN THE ENDS OF THE 



STACK. 

toz= (tempR (1, k) -tempL (1, k) ) /lengsec (1, k) ; % approx 

dz/dt 

ns=0; 

for la=k+numblay (1 , k) -1 : -1 : k 
ns=ns+l ; 

tens=tempR ( 1 , k) - (tempR ( 1 , k) -tempL ( 1, k) ) *ns/numblay ( 1 , k) ; 
tnsml=tempR ( 1 , k) - (tempR ( 1 , k) -tempL (l,k) ) * (ns-1) /numb lay ( 1 , k) 

r 

tave (1, la) = (tens+tnsml) / 2; 

dsub (1, k) =lengsec (1, k) /numblay (1, k) ; 

dens (1, la) =ampres*4 . Oe-3 / (tave (1, la) *8. 3143) / % density 

visc(l,la)=1.887e-5* (tave (l,la)/273.15) / '(.6567); % 

viscosity 

kgas (1, la) =visc (1 , jup) *cp/npr; % termal conductivity 
sspeed (1, la) =972 . 8*sqrt (tave (1, la) /273 . 15) ; % speed 

% isothermal 

% GET LAMBDA AND LAMBDA (T) 

lambda (l,la)=ratio(l,k) . *sqrt (dens (1,1a) . *w. /vise (1, la) ) ; 
lambdt (1, la) =sqrt (npr) * lambda (1, la) ; 
ar ( 1 , la) =sqrmi * lambda (1,1a) ; 
argum (1, la) =exp (-ar (1, la) ) / 
ar (1, la) =ar (1, la) /2; 

ctanh ( 1 , la) = ( 1 -argum (1,1a)) . / ( 1+argum (1,1a)); 
flam (1, la) =l-ctanh (1, la) . /ar (1, la) ; % f ( 1 ) 



59 



ar (1, la) =sqrmi*lambdt (1, la) ; 
argum (1, la) =exp (-ar (1, la) ) ; 
ar (1, la) -ar (1, la) 12 ; 

ctanh (1, la) = (1-argum (1,1a))./ (1 + argum (1, la) ) ; 
f lamt (1, la) =l-ctanh (1, la) . /ar (1, la) ; % f(l(t)) 

low=w/sspeed (1, la) ; 

ka (1, la) =low. *sqrt ( (gamma- (gamma-1) . . . 

. *f lamt (1,1a) ) ./flam (1,1a) ) ; 
zint (1,1a) =dens (l,la).*w... 

. / (poros ( 1 , k) . *f lam ( 1 , la) . *ka (1 , la) ) ; 
alprim (1, la) =toz* (flamt (1, la) . . . 

. /flam (1, la) -1) . / (2*tave (1, la) * (1-npr) ) ; 
kal=ka (1, la ) . ' ; 
zintl=zint (1, la) . ' ; 
alprl=alprim (1, la) . ' ; 

zetaO= (lengsec (1, k) - (ns-1) *lengsec (1, k) /numblay (1, k) ) ; 
zetafin= (lengsec (1, k) - (ns) * lengsec (1, k) /numblay (1, k) ) ; 
zO=z (1, 1 + la ) . ' ; 
tol=l . e-4 ; 

[ zeta, zout ] =trode45 ( ' tube' ,zetaO,zetafin,zO,tol,kal, zint 1, al 
prl) ; 

n=length (zout) ; 
z ( 1 , la) =zout (n, 1 ) . ' ; 
end 
end 
end 
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z=z (1,1); 
dens=dens (1,1); 
sspeed=sspeed (1,1); 
visc=visc (1,1); 
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